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We present Monte Carlo simulations on a new class of lattice models in which the degrees of 
freedom are elements of an abelian or non-abelian finite symmetry group Q, placed on directed 
edges of a two-dimensional lattice. The plaquette group product is constrained to be the group 
identity. In contrast to discrete gauge models (but similar to past work on height models) only 
elements of symmetry-related subsets S G Q are allowed on edges. These models have topological 
sectors labeled by group products along topologically non-trivial loops. Measurement of relative 
sector probabilities and the distribution of distance between defect pairs are done to characterize 
the types of order (topological or quasi-LRO) exhibited by these models. We present particular 
models in which fully local non-abelian constraints lead to global topological liquid properties. 



Introduction : Ever since Wen first proposed that 
"Topological Order" be considered as a means of classify- 
ing chiral spin states in superconductors pQ, applications 
of this "post-Landau" paradigm have caused great excite- 
ment in applications to quantum systems — the quantum 
Hall effect [2], fractional quantum Hall effect [3J, spin 
liquids [IHZ] , and topological quantum computation[5]. 
Most other physical phenomena, such as critical phe- 
nomena, were first understood classically prior to being 
realized in quantum mechanics, so the question arises 
whether there is a useful notion of topological order in a 
classical setting. Let us define a generalized topological 
order in a classical ensemble by the existence of sectors 
of (energetically or entropically) degenerate states, dis- 
connected in the thermodynamic limit, which cannot be 
distinguished by any local order parameter. Indeed, one 
can realize this sort of order in classical models where 
the inaccessibility is built in by hand, or in a limit where 
inter-sector transitions are suppressed by an activation 
energy much larger than the temperature [9l [TUj. The 
motivation for such models is that first, by sharing many 
features with the quantum models, they provide an arena 
to do calculations that would not be feasible in the quan- 
tum case; secondly, such a classical ensemble can furnish 
the Hilbert space for contructing a quantum model with 
the same topological order (e.g. the classical Z 2 topolog- 
ical order in dimer coverings on triangular lattices |11|). 

As a particular instance of classical topological order- 
ing, one of us has proposed |10j a new family of models 
which we will henceforth call "generalized height mod- 
els" which are based on an abelian or non-abelian dis- 
crete group in the same way that lattice "height" models 
[T21 - H6] are based on the integers. These models share 
many properties of their quantum analogues, including 
massively degenerate ground states, topological defect 
charges, and the possibility in non-abelian cases of com- 
bining two defect charges in more than one way (called 
"fusion channels" in the quantum context). In this 
paper we present simulation results from these models, 
focused primarily on two measurements which charac- 



terize whether a given model is topologically ordered: 
the relative probabilities of being each sector, and the 
distribution of separations between a pair of topological 
defects. 

Simulations of the Model: An extended introduction 
of generalized height models was given previously [10] : 
here we will outline the properties of a general model 
and discuss the specific models we have already looked 
at. Generalized height models can be thought of as lying 
between height models [I2HT6] and discrete lattice-gauge 
models |17j . The former are models in which the vertices 
of a (not necessarily Bravais) lattice are labeled by finite 
integer heights (or, conversely, directed edges are labeled 
by the integer differences between those heights) ; whereas 
the latter are defined by elements of general groups as- 
signed to the edges of a lattice (we will later call them 
"spins"). 

We can define them via two constraints. These mod- 
els are generalizations of height models in that the set of 
integers used in height models to label the degrees of free- 
dom on the edges are replaced by any discrete group (fi- 
nite or infinite) Q as in gauge models. The "well-defined" 
constraint of integers around a plaquette adding to zero 
is thus replaced with the constraint on edge elements 
<r(ri,iv) on edges between vertices at and rv which 
defines our massively degenerate ensemble: 

JJ a(T U r i+1 ) = eg. (CI) 

This constraint is wholly local, and can be abelian or non- 
abelian given the group Q. If this constraint is violated 
on a plaquette, that plaquette is said to have a charge 
defect, with the charge being the plaquette product. 

The second constraint is the true novelty of general- 
ized height models. Unlike discrete gauge models, but 
like height models, only elements of a chosen symmetry- 
related subset 

Scg (C2) 



are placed on the directed edges of the lattice. The set 
S is normally a collection of symmetry classes, and is 
strictly not equal to Q. Essentially, S is a symmetry- 
related subset of elements in Q. 

In discrete gauge models, all group elements are al- 
lowed and therefore spins living on lattice edges may not 
be correlated with each other (beyond the gauge con- 
straint) as in the toric code [TS], for example. The cor- 
relations between spins in these restricted spin subset 
models may be exponentially decaying, but are generally 
non-trivial. In fact, by changing the size of the allowed 
spin subset we can discretely interpolate between totally 
free gauge models and quasi-long range ordered states 



Note that constraint CI implies that elements along 
any topologically trivial loop in the lattice must multiply 
to the identity, but it does not imply that a topologi- 
cally non-trivial loop need do so. Our simulations take 
place on the torus, and we therefore have two indepen- 
dent non-trivial loop products which we take as labels for 
the disjoint partitions of our ensemble, hereafter referred 
to as sectors. In the case of abelian Q, each sector has a 
uniquely defined label ; for non-abelian groups they are 
defined up to conjugacy[TD]. Local updates (defined 
precisely below) cannot take us from a configuration in 
one sector to a configuration in another, but cntropically 
expensive global updates can. Monte Carlo simulations 
involving updates of these two types provide us with a 
means of measuring the entropic cost of each sector. 
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1. A sample configuration of our model with group 
'2 x Tj-z and spin sub-set S = Z2 x Z2 \ e on a periodic 
square lattice. The upper left plaquette is a "normal" plaque- 
tte: the product of the edges is the identity. The lower left and 
upper middle plaquettes (connected via the black dashed line) 
are defect plaquettes with opposite "charges" (non-identity 
plaquette products). We have also shown the plaqutte prod- 
uct along one of the topologically non-trivial loops — this will 
be half of the "sector" label for this configuration. 

In our Monte Carlo simulations, a sequence of local and 
global update moves (satisfying detailed balance) is used 
to transition between different configurations. The local 
move is a single-vertex update in which we multiply all 



outgoing edges for one vertex by some group value g € Q. 
This update preserves the plaquette and sector products 
(up to conjugacy), but is not always allowed (if any of 
the resulting elements on the outgoing edges are not in 
the spin set <S, that update is rejected). 

Transitioning between sectors requires a global move. 
Our chosen move involves creating a defect / anti-defect 
pair, randomly walking one around the lattice, and let- 
ting them fuse into a new defect (non-abelian models 
only) or annihilate. The entire process is referred to as a 
completed defect walk, and an individual jump of a defect 
from one plaquette to the neighbor as a defect step. 

In the case of an Abelian symmetry group Q the 
defect / anti-defect pair will always annihilate upon 
coming back together, but in models generated from 
Non- Abelian symmetry groups they can potentially re- 
combine into a single defect of a different charge. If a 
defect walks once around some topologically non-trivial 
loop, every transverse loop gets crossed once and has its 
loop product changed, so that the new configuration re- 
sides in a different sector. Such updates, that change sec- 
tors while satisfying detailed balance, allow us to measure 
the relative weights of the different sectors and thereby 
ascertain whether the ensemble is topologically ordered. 

Constraint IC2I induces non-trivial correlations be- 
tween spins in the lattice, and therefore walking defects 
may interact with each other through the intermediary 
spins. Exponentially decaying probability distributions 
of the distance between defect pairs indicate topologi- 
cally ordered liquid-like phases with exponentially decay- 
ing correlations and deconfined topological defects. Even 
though we adopt the ensemble of maximum entropy, it 
is possible that the resulting ensemble has more than 
topological order [TO]: e.g. constraints (CI) and (C2) 



together may so constrain the possibilities in each pla- 
quette that we get a height model, or we could find an 
emergent long-range ordered phase. 

We focus primarily on the square lattice with groups 
and spin subsets as enumerated in table |TJ To gauge the 
difference between models with topological order or with 
more order, we decided to use the model Q = Z5 with 
S = {±1} — which is in a one-to-one mapping with 
the well-studied six-vertex model (also known as the ice 
model [5D]) — as a benchmark. We ensured that sim- 
ulations on this model reproduced known values for the 
six-vertex model before proceeding with possibly topo- 
logically ordered models. The other group choices were 
motivated for different reasons: a model using the three 
elements of the abelian group Z 2 x Z 2 can behave like 
a three-coloring model in some situations; the group S3 
is the smallest non-abelian group; and A5 is both the 
smallest non-abelian simple group (meaning it cannot be 
"factored" into quotient groups) and is also the symme- 
try group of the regular icosahedron. 
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Group Q 


Spin set S 


A or Non-A 


Z 5 


{+1,-1} 


A 


Z 2 x Z 2 


{(0,1), (1,0), (1,1)} 


A 




All non-identity elements 


NA 


A 


All elements of order three 


NA 



TABLE I. Models simulated for this paper. 



In summary, generalized height models offer us two 
different "tunable" parameters: the local symmetry of 
the model, as defined by the full group Q ; and the subset 
of allowed spins. Each different choice of spin set S can 
be thought of as tuning a family of models [19] . 

Our measurements are based on Monte Carlo simula- 
tions of the models with group choices as described in 
table |T] and periodic lattices of edge length L. Define a 
"sweep" of single-site updates as L 2 consecutive single- 
vertex updates at randomly selected vertices. Global de- 
fect walks are performed one step at a time, with the ratio 
of single-site sweeps to individual steps in the defect ran- 
dom walk specified by hand (for most of the simulations 
presented in this paper that ratio was 1 sweep per 100 
defect steps). We have performed simulations on lattice 
sizes ranging from 4 2 to 256 2 sites. 

Sector Probabilities : By definition, the primary prop- 
erty of topological order in our models is an entropic 
degeneracy of the lowest energy ensemble in the ther- 
modynamic limit. To demonstrate this, we need some 
measurement of the relative number of states in different 
sectors. Completed defect walks traversing a topologi- 
cal^ non-trivial path through the system will transition 
our configuration into a different sector. The relative 
probability P(T) that we would find ourselves in sector 
r = (7,7') after such a walk then is proportional to the 
(relative) number of configurations in a sector. A topo- 
logically ordered state will have P(T) converging to a 
constant value P^ independent of T as the lattice size 
L increases; a critical state will have P(T) converging to 
different values for different V. 

In contrast, for models with topological order we ex- 
pect the probability of being in a sector T at any lattice 
length L to depend on the class, and for those probabil- 
ities to decay to each other as a function of lattice size, 



P(T,L) = P oo (T) + b L P (T)e^(-L/£(T)), (1) 

where b L is a prefactor term which only depends on the 
even- or odd-ness of L (often b L = (— 1) L ). By defini- 
tion, for a model to be topologically ordered, the number 
of configurations in each of the N s sectors must be the 
same in the thermodynamic limit: Poo(r) = 1/N S for all 

r. 

We start with the benchmark six-vertex model: the 
number of configurations per sector T = (7, V)[2T] in 



the six-vertex model can be calculated analytically 
P(r,i)cxexp(-|( 7 2 +7' 2 )), « 
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(2) 



Note that for any integer height model such as this 
P(T,L) is not dependent on the lattice edge size L [22] . 
The constancy as a function of L can be seen in figure 
[2] and the dependence on the sector T in the subset 
of that figure. Our simulations reproduced a value of 
k = . 523(1). 




FIG. 2. For the model S3 (2, 3) (solid lines), the number 
of configurations in a sector converges to the same constant 
value as a function of lattice size L. Though all sectors con- 
verge to the same number, at short lattice lengths they behave 
differently up to the type of group element (i.e. the orders of 
the elements labeled in the figure). For Z5 (dashed lines), the 
probability converges to a constant which depends on the sec- 
tor T as P(T,L) oc exp(-(K/2)(7 2 + 7 ' 2 )). This exponential 
scaling of sector probabilities for Z5 is shown in the inset. 



The foreground of figure [2] shows the sector probability 
of the non-abelian £3(2, 3)sg model. The sector probabil- 
ities for different sectors T converge to a shared constant 
value, indicating topological order. The data plotted in 
figure[2]was averaged over different sectors that are equiv- 
alent by symmetry and thus the curves there are labeled 
not by the group elements of particular sectors but rather 
by the symmetry class of the element (e.g. "2" for order- 
two elements, "e" for the identity, etc.). 

Fit values for different sectors for some sectors for a 
few select models are compiled in table [TTJ Caution: sim- 
ulations of the A5 model are not fully understood yet, 
though sector probabilities do appear to be decaying to 
the same constant value. 



Defect Pair Distributions: The constraint ( C2 1 gener 



ates non-trivial local correlations. We expect to see ex- 
ponentially decaying correlations in the topologically or- 
dered models and power law correlations in the Z5 model. 
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-Po(r) 


£(T) 


r 


W) 


*(r) 


Q = 


Z 2 x Z 2 ,S 


= G\e 


Q = A^,S = Order 2 elem. 


(e, e) 


.3138(8) 


1.9290(2) 


(p p) 


.184(5) 


1.63(2) 


(a e) 


.0431(1) 


2.749(4) 


(p 2) 


0233(6"l 


1.96(2) 


(a a) 


.01679(9) 


4.45(2) 


(p 3") 


01 09f?^ 




(a, a') 


.1209(4) 


1.835(3) 


(e, 5i) 


.0022(8) 


3.0(5) 


g 




g\ e 


(e,5 2 ) 


.002(1) 


3.1(9) 


(e, e) 


.45(2) 


0.838(7) 


(2,2) 


.048(3) 


1.06(2) 


(2,e) 


.061(7) 


0.95(3) 


(3,3) 


.0122(5) 


1.49(2) 


(3,e) 


.053(3) 


0.96(1) 


(5i,5i) 


.0131(3) 


2.20(3) 


(2,2) 


.097(6) 


0.98(1) 


(5 2 ,5 2 ) 


.0131(3) 


2.20(3) 


(3,3) 


.20(2) 


0.86(1) 


(5i,5 2 ) 


.0138(2) 


2.21(2) 



TABLE II. Fit parameters for the assumed exponential fit- 
ting form |l]) parameters for Abelian and non-Abelian models. 
The labels F are grouped by conjugacy class (not individual 
element) as sectors are only defined modulo the group conju- 
gacy. The sectors for the S3 and ,4s simulations are labeled 
by the order of the conjugacy class (^5 conjugacy classes 5i 
and 52 are related by an outer automorphism), whereas the 
labels for the Z2 x Z2 model are a for any of the three non- 
identity elements of that group, and (a, a') when 7 and 7' are 
not the same element. 




FIG. 3. The left panel shows the defect pair distribution 
function -F(r) for the only defect / anti-defect type possible 
in the Zs{±l}sq model (+1 and —1) for many different lat- 
tice sizes L. The plot is on a log-log scale, clearly indicating 
the power law behavior. The right panel depicts the defect 
pair distribution function for the defect types in the Z2 x Z2 
model at only L = 24. The exponential decay is manifested 
in a rapid decay to a constant background value, indicating 
deconfinement of defects. 



For our topologically ordered models, the defects are de- 
confined at large distances, but have considerable inter- 
actions at short ranges. 

The distribution function Fgj' (r) of the distances be- 
tween any two defects S, 5' defines an effective entropic 
potential V(r) via 

F StS ,(T) = F Q exp(-V s , s >(r)). (3) 

For the six- vertex model (here Z5), we may only 
have defects of charge ±1, and furthermore it is known 
that these defects behave like Coulomb charges in two 
dimensions [TJJ [23] • The energy required to hold them 
separate is logarithmic 

V(r) -(2/7r)«Iog(r), (4) 

with k = 7r/6 as noted before , so the pair distribution 
follows a power law: 

F ±2 , ±2 (r)cx±exp(^log(r)^ = r^ 3 . (5) 

The power-law behavior can be seen in the distribution 
functions for the Z5 model for several lattice sizes L 
shown in the left panel of figure [3] 

For topological liquids, on the other hand, exponen- 
tially decaying correlations result in a potential which 
decays exponentially to a constant (that is to say, defects 
are deconfined) . This behavior for the Z 2 x Z 2 model is 
shown on the right side of figure [3j 

Analytical transfer matrix toy calculations [TU] suggest 
that the form of pi) ought to be a sum of exponentially 



decaying terms. Therefore, we fitted our results to the 
functional form 

log (F StS , (r)) = F%, + exp (-r/&,y) , (6) 

with the fit parameters being given in Table |III| for the 
models which show the exponentially deconfined defect / 
anti-defect pair distributions characteristic of TO mod- 
els. 
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s 


S, 6' 




%S,8' 


Z 2 x Z 2 


Z 2 x Z 2 \ e 


(a, a) 


0.170(6) 


0.54(2) 


s 3 


S 3 \e 


(2,2) 
(3,3) 


0.33(9) 
0.11(5) 


0.20(1) 
0.26(4) 



TABLE III. TO models feature exponentially decaying defect 
pair distributions between defect types S and 5' . Fits were 
performed on the TO models simulated for this paper and 
the functional form that was used for the fitting is given in 
equation (|6| . For the model Z2 x Z2 we have averaged all three 
defect / anti-defect pairs; for S3 we have averaged all pairs of 
order 2 and order 3 defect / anti-defect pairs, respectively. 



Conclusions : Using a generalized definition of topo- 
logical order, we have explored a new family of classical 
models which exhibit properties analogous to quantum 
systems with quantum topological order: (1) topologi- 
cally robust partitions of massively degenerate ensem- 
bles, (2) finite size effects, (3) defect deconfinement (in 
topological systems), and (4) mediated interactions be- 
tween charged defects. Elsewhere, we will also exhibit (5) 
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an ability to "tune" the number of degrees of freedom 
to transition between topologically ordered and quasi- 
critical states [15]. 

Though the applications of topological order are be- 
coming more and more widespread, the number of 
tractable models that realize topological order is un- 
fortunately still limited. Our generalized height mod- 
els expand this collection. Furthermore, this non- 
perturbative approach may be illuminating in unrav- 
eling the phenomenology of condensed matter systems 
where a gauge symmetry is realized only in perturba- 
tive limits. Various sorts of quantum models could be 
constructed from these models (either via the Rokhsar- 
Kivelson construction 24-26] or the Levin- Wen "string- 
net" construction [57] Whether the non-abelian models 
contain anyonic defects depends on the phase factors as- 
signed to the matrix elements in quantizing the general- 
ized height model in question. 
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